### GENERATING FIGURE C.1: LIST EXPERIMENTS CONDITIONAL ON ASPIRATION-ABILITY GAP

### USES PNAS DATASET

### REQUIRED PACKAGES
require(multiwayvcov)
require(lmtest)
require(dplyr)
require(gridExtra)
require(grid)
require(ggplot2)


#running conditional models with pnas
work.gap.aspab <- lm(worklist.pre ~ worklist.pre.dum + aspab.gap + worklist.pre.dum*aspab.gap 
                     + gender + age + education + income.imp + freq.help.net + factor(parish), 
                     data= pnas)

summary(work.gap.aspab)


##generating clustered standard errors
work.gap.aspab.clust<-cluster.vcov(work.gap.aspab, pnas$ed.parish)
work.gap.aspab.mcse <- coeftest(work.gap.aspab, work.gap.aspab.clust)
work.gap.aspab.mcse


## Pull out data for marginal effects plots
# coefficients
work.gap.aspab.b2 <- work.gap.aspab.mcse[2,1]
work.gap.aspab.b3 <- work.gap.aspab.mcse[3,1]
work.gap.aspab.b22 <- work.gap.aspab.mcse[22,1] # interaction term
# variances
work.gap.aspab.varb2 <- work.gap.aspab.clust[2,2]
work.gap.aspab.varb3 <- work.gap.aspab.clust[3,3]
work.gap.aspab.varb22 <- work.gap.aspab.clust[22,22]
# covariances
work.gap.aspab.covb2b22 <- work.gap.aspab.clust[2,22]
work.gap.aspab.covb3b22 <- work.gap.aspab.clust[3,22]

list(work.gap.aspab.b2, work.gap.aspab.b3, work.gap.aspab.b22, work.gap.aspab.varb2, 
     work.gap.aspab.varb3, work.gap.aspab.varb22, work.gap.aspab.covb2b22, work.gap.aspab.covb3b22)

##transforming dataframe to graph
work.aspab.gap.df <- pnas %>% mutate(aspab.gap = seq(0, 6, length.out=1163),
                                        work.gap.aspab.conbx = work.gap.aspab.b2+work.gap.aspab.b22*aspab.gap,
                                        work.gap.aspab.consx = sqrt(work.gap.aspab.varb2 + work.gap.aspab.varb22*(aspab.gap^2)+
                                                                      2*work.gap.aspab.covb2b22*aspab.gap),
                                        work.gap.aspab.ax = 1.96*work.gap.aspab.consx,
                                        work.gap.aspab.upperx = work.gap.aspab.conbx + work.gap.aspab.ax,
                                        work.gap.aspab.lowerx = work.gap.aspab.conbx - work.gap.aspab.ax)

##interaction plot with confidence intervals based on clustered SEs
work.aspab.gap.plot <- ggplot(data = work.aspab.gap.df) + 
  scale_y_continuous("Estimated Support for Fake Documents") + 
  scale_x_continuous(limits = c(0,6)) +
  geom_line(aes(x = aspab.gap, y = work.gap.aspab.conbx)) + 
  geom_ribbon(aes(x=aspab.gap, ymin=work.gap.aspab.lowerx, ymax = work.gap.aspab.upperx), fill="red", alpha=.2) + 
  geom_point(aes(x=aspab.gap, y=work.gap.aspab.conbx), size = 1, color = "red") + theme(legend.position="none") +  
  geom_hline(yintercept = 0) +
  xlab("Ability gap") +
  ggtitle("Work Under Table") + theme(plot.title = element_text(face="bold"))
work.aspab.gap.plot


enter.gap.aspab <- lm(enterlist.pre ~ enterlist.pre.dum + aspab.gap + enterlist.pre.dum*aspab.gap 
                      + gender + age + education + income.imp + freq.help.net + factor(parish), 
                      data= pnas)

summary(enter.gap.aspab)

##generating clustered standard errors
enter.gap.aspab.clust<-cluster.vcov(enter.gap.aspab, pnas$ed.parish)
enter.gap.aspab.mcse <- coeftest(enter.gap.aspab, enter.gap.aspab.clust)
enter.gap.aspab.mcse
rownames(enter.gap.aspab.mcse)[2] <- "Treat"
rownames(enter.gap.aspab.mcse)[22] <- "Treat*Aspiration Ability Gap"



## Pull out data for marginal effects plots
# coefficients
enter.gap.aspab.b2 <- enter.gap.aspab.mcse[2,1]
enter.gap.aspab.b3 <- enter.gap.aspab.mcse[3,1]
enter.gap.aspab.b22 <- enter.gap.aspab.mcse[22,1] # interaction term
# variances
enter.gap.aspab.varb2 <- enter.gap.aspab.clust[2,2]
enter.gap.aspab.varb3 <- enter.gap.aspab.clust[3,3]
enter.gap.aspab.varb22 <- enter.gap.aspab.clust[22,22]
# covariances
enter.gap.aspab.covb2b22 <- enter.gap.aspab.clust[2,22]
enter.gap.aspab.covb3b22 <- enter.gap.aspab.clust[3,22]

list(enter.gap.aspab.b2, enter.gap.aspab.b3, enter.gap.aspab.b22, enter.gap.aspab.varb2, 
     enter.gap.aspab.varb3, enter.gap.aspab.varb22, enter.gap.aspab.covb2b22, enter.gap.aspab.covb3b22)

##transforming dataframe to graph
enter.aspab.gap.df <- pnas %>% mutate(aspab.gap = seq(0, 6, length.out=1163),
                                         enter.gap.aspab.conbx = enter.gap.aspab.b2+enter.gap.aspab.b22*aspab.gap,
                                         enter.gap.aspab.consx = sqrt(enter.gap.aspab.varb2 + enter.gap.aspab.varb22*(aspab.gap^2)+
                                                                        2*enter.gap.aspab.covb2b22*aspab.gap),
                                         enter.gap.aspab.ax = 1.96*enter.gap.aspab.consx,
                                         enter.gap.aspab.upperx = enter.gap.aspab.conbx + enter.gap.aspab.ax,
                                         enter.gap.aspab.lowerx = enter.gap.aspab.conbx - enter.gap.aspab.ax)

##interaction plot with confidence intervals based on clustered SEs
enter.aspab.gap.plot <- ggplot(data = enter.aspab.gap.df) + 
  scale_y_continuous("Estimated Support for Fake Documents") + 
  scale_x_continuous(limits = c(0,6)) +
  geom_line(aes(x = aspab.gap, y = enter.gap.aspab.conbx)) + 
  geom_ribbon(aes(x=aspab.gap, ymin=enter.gap.aspab.lowerx, ymax = enter.gap.aspab.upperx), fill="red", alpha=.2) + 
  geom_point(aes(x=aspab.gap, y=enter.gap.aspab.conbx), size = 1, color = "red") + theme(legend.position="none") +  
  geom_hline(yintercept = 0) +
  xlab("Ability gap") +
  ggtitle("Fake Documents") + theme(plot.title = element_text(face="bold"))
enter.aspab.gap.plot

#enter and work aspab gap graph
aspab.gap.df <- cbind(work.aspab.gap.df, enter.aspab.gap.df)
aspab.gap.plot.combined <- ggplot(data = aspab.gap.df) + 
  theme(text = element_text(size=15),
        plot.title = element_text(size =10)) +
  scale_y_continuous("Estimated Support") + 
  scale_x_continuous(limits = c(0,6)) +
  geom_ribbon(aes(x=aspab.gap, ymin=work.gap.aspab.lowerx, ymax = work.gap.aspab.upperx), fill="blue", alpha=.2) + 
  geom_line(aes(x=aspab.gap, y=work.gap.aspab.conbx, color = "Work Undertable"), size = 1) +  
  geom_ribbon(aes(x=aspab.gap, ymin=enter.gap.aspab.lowerx, ymax = enter.gap.aspab.upperx), fill="red", alpha=.2) + 
  geom_line(aes(x=aspab.gap, y=enter.gap.aspab.conbx, color = "Fake Documents"), size = 1)  +  
  geom_hline(yintercept = 0) +
  xlab("Ability Gap") + theme(legend.position="bottom") +
  ggtitle("Ability Gap and Support for Illegal Strategies") + theme(plot.title = element_text(face="bold")) +
  scale_colour_manual(name='',values=c('Work Undertable'='blue','Fake Documents'='red'))
aspab.gap.plot.combined
